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Abstract 

We study the kinematics of multigrid Monte Carlo algorithms by 
means of acceptance rates for nonlocal Metropolis update proposals. An 
approximation formula for acceptance rates is derived. We present a 
comparison of different coarse-to-fine interpolation schemes in free field 
theory, where the formula is exact. The predictions of the approxima- 
tion formula for several interacting models are well confirmed by Monte 
Carlo simulations. The following rule is found: For a critical model with 
fundamental Hamiltonian 7i(0), absence of critical slowing down can 
only be expected if the expansion of (7i(0 + ^)) in terms of the shift ip 
contains no relevant (mass) term. We also introduce a multigrid update 
procedure for nonabelian lattice gauge theory and study the acceptance 
rates for gauge group SU (2) in four dimensions. 
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1 Introduction 



Monte Carlo simulations of critical or nearly critical statistical mechanical systems 
with local algorithms suffer from critical slowing down (CSD). Roughly speaking, 
the autocorrelation time in the Markov chain behaves like r ~ in the vicinity 
of a critical point, where ^ denotes the correlation length, and z is the dynamical 
critical exponent. For conventional local algorithms, z 1. For accelerated local 
algorithms such as overrelaxation or the optimized hybrid Monte Carlo algorithm, 
one can sometimes achieve ^ ~ 1 [|T], 0. To overcome the problem of CSD, various 
nonlocal Monte Carlo algorithms have been developed. 

Cluster algorithms 0] are successful in overcoming CSD for a large class of mod- 
els. The alternative is multigrid Monte Carlo @, |§, ^ . In this paper, every algorithm 
that updates stochastic variables on a hierarchy of length scales is called multigrid 
Monte Carlo algorithm. There are models where no successful cluster algorithms 
have been found whereas multigrid Monte Carlo algorithms work ||^, 

Presently, the only generally applicable method to study algorithms for interact- 
ing models is numerical experiment. It is however important to have some theoretical 
understanding that helps to predict which algorithms will have a chance to over- 
come CSD in simulations of a given model. As a contribution to the research in 
this direction we here present a study of the kinematics of multigrid Monte Carlo 
algorithms With kinematics we here mean the study of the scale (block size) 
dependence of the Metropolis acceptance rates for nonlocal update proposals. We 
do not address the much more ambitious problem of analytically investigating the 
full dynamical critical behavior of the stochastic processes involved. Our analysis is 
nonetheless of relevance because sufficiently high acceptance rates are necessary for 
multigrid Monte Carlo procedures to overcome CSD. 

We derive an approximation formula for the block size dependence of acceptance 
rates for nonlocal Metropolis updates. The influence of the coarse-to-fine interpola- 
tion kernel (shape function) on the kinematics in free field theory, where the formula 
is exact, is investigated in detail. 

The formula is then applied on several interacting models and turns out to be 
a very good approximation. We find necessary criteria for a given multigrid algo- 
rithm to eliminate CSD: For a critical model with a fundamental Hamiltonian 7^(0) 
absence of CSD can only be expected if the expansion of (7^(0 -t- '0)) in terms of the 
shift contains no relevant term (mass term). 

There is an urgent demand for accelerated Monte Carlo algorithms in lattice 
gauge theory. The present state-of-the-art algorithm is overrelaxation How- 



ever, effort was also spent in developing nonlocal algorithms for gauge theories. An 
efficient cluster algorithm was found for SU (2) gauge theory at finite temperature, 
however only in the special case A^j = 1 [O. For a recent cluster algorithm ap- 



proach to f/(l) gauge theory see [Q. Multigrid algorithms for t/(l) gauge models 



were introduced and studied in two and four dimensions [1^, |14|. A different but 



^ Parts of this paper are published in short form in 
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related nonlocal updating scheme in the abelian case is the multiscale method |T5|. 

In this paper, we propose a multigrid algorithm for nonabelian gauge theory and 
analyze its kinematics. Our approximation formula turns out to be very reliable 
also in this case and allows for a prediction of acceptance rates for a large class of 
nonlocal updates. 

We believe that the proposed algorithm should be able to accelerate the local 
Monte Carlo dynamics at least by a constant factor. 

This paper is organized as follows: In section ^ we introduce multigrid Monte 
Carlo algorithms. Section |^ contains the derivation of our approximation formula for 
acceptance rates. Several coarse-to-fine-interpolation kernels are discussed in section 
^. In section ^ the acceptance rates in free field theory are examined in detail. The 
kinematical analysis for the Sine Gordon, XY and 0"^ models is presented in section 
^. In section |^ we propose a multigrid procedure for nonabelian gauge theories and 
analyze its kinematics. A summary is given in section & 



2 Multigrid Monte Carlo algorithms 

We consider lattice models with partition functions 

Z = j X{dct>x exp(-7^(0)) (1) 

on cubic (i-dimensional lattices Aq with periodic boundary conditions. The lattice 
spacing is set to one. We use dimensionless spin variables An example is single- 
component 0^-theory, defined by the Hamiltonian 

n{<p) = 1(0, -A0) + ^ E 0^ + ^ E 0^ (2) 

where 

(0,-A0)= E (0.-0.)'- (3) 

<x,y> 

The sum in eq. is over all nearest neighbor pairs in the lattice 0. 

A standard algorithm to perform Monte Carlo simulations of a model of the 
type defined above is the local Metropolis algorithm: One visits in a regular or 
random order the sites of the lattice and performs the following steps: At site Xo, 
one proposes a shift 

K ^ 0L„ = K + s. (4) 

The configuration {0^;} remains unchanged for x ^ Xq. 5 is cl random number 
selected according to an a priori distribution p{s) which is symmetric with respect 
to s — > —s. E.g., one selects s with uniform probability from an interval [— e,^]. 
One then computes the change of the Hamiltonian 

AH = H{(j)') - H{(j)) . (5) 

^The definitions for lattice gauge theory will be introduced in section M 
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Finally the proposed shift is accepted with probability min[l, exp(— A7i)]. Then one 
proceeds to the next site. 

The local Metropolis algorithm suffers from CSD when the correlation length 
in the system becomes large: long wavelength fluctuations cannot efficiently be 
generated by a sequence of local operations. It is therefore natural to study nonlocal 
generalizations of the update procedure defined above. 

Consider the fundamental lattice Aq as divided in cubic blocks of size l'^. This 
defines a block lattice Ai. By iterating this procedure one obtains a whole hierarchy 
of block lattices Aq, Ai, . . . , A^^ with increasing lattice spacing. This hierarchy of 
lattices is called multigrid. 

Let us denote block lattice points in A^ by x'. Block spins are defined on 
block lattices A^. They are averages of the fundamental field (px over blocks of side 
length Lb = I'': 

^x' = Lt''^"L~^'Y.<P^- (6) 

Xdx' 

The sum is over all points x in the block x' . The L^-dependent factor in front of 
the average comes from the fact that the corresponding dimensionful block spins are 
measured in units of the block lattice spacing: A scalar field </'(x) in d dimensions 
has canonical dimension (2 — d)/2. Thus = a^^"'^^/^^^., where a denotes the 
fundamental lattice spacing. Now measure the dimensionful block spin in units 
of the block lattice spacing a': = a'(^~'^)/^$j./ , with a' = aLs- If we average 

in a natural way = L'^'^ Yl,x&x' 4>{^) aiid return to dimensionless variables, we 

obtain eq. (|]). 

A nonlocal change of the configuration (p consists of a shift 

(l)x^ (l)x + S^Jx- (7) 

s is a real parameter, and the "coarse-to-fine interpolation kernel" (or shape func- 
tion) ipx determines the shape of the nonlocal change, is normalized according 
to 

L^'^ Y.i,x = Lt'^''5x'^x'^ . (8) 

Note that by the nonlocal change (|^, the block spin is moved as $2,/ $2,/ + s for 
x' = Xg, and remains unchanged on the other blocks. The simplest choice of the 
kernel ijj that obeys the constraint (P) is a piecewise constant kernel: ipx = L^b'^^^"^^ 
if x G x'q, and otherwise. Other kernels are smooth and thus avoid large energy 
costs from the block boundaries. A systematic study of different kernels will be 
given in section ^ below. 

The s-dependent Metropolis acceptance rate for such proposals is given by 

n{s) = (min[l,exp(-A7^)]) . (9) 

Here, ((.)) denotes the expectation value in the system defined by eq. (|l|). Further- 
more, 

AH = ni(i) + stp) - n{(i)) . (10) 
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fl{s) can be interpreted as the acceptance rate for shifting block spins by an amount 
of s, averaged over a sequence of configurations generated by a Monte Carlo simu- 
lation. Note that fl{s) does not depend on the algorithm that we use to compute 
it. Q{s) is a useful quantity when one wants to know how efficiently updates with 
increasing nonlocality (i.e. increasing block size Lb) can be performed. Of course, 
different choices of the kernel ijj result in different acceptance rates. 

In actual Monte Carlo simulations, s is not fixed. In the same way as in the 
local Metropolis algorithm, s is a random number distributed according to some a 
priori probability density. If we choose s to be uniformly distributed on the interval 
[—e, e], the integrated acceptance rate Pace (as customarily measured in Monte Carlo 
simulations) is obtained by averaging Q{s) as follows: 



It turns out to be a good good rule to adjust the maximum Metropolis step size e 
such that Pacci^) ~ 50%. 

We consider every algorithm that updates stochastic variables on a hierarchy of 
length scales as multigrid Monte Carlo algorithm. However, there are two different 
classes of multigrid algorithms: multigrid algorithms in a unigrid implementation 
and "explicit" multigrid algorithms. 

In the unigrid formulation one considers nonlocal updates of the form (0). Up- 
dates on the various layers of the multigrid are formulated on the level of the finest 
lattice Aq. There is no explicit reference to block spin variables $ defined on coarser 
layers with k > 0. In addition, unigrid also refers to a computational scheme: 
Nonlocal updates are performed directly on the level of the finest grid Aq in practical 
simulations. 

In contrast, the explicit multigrid formulation consists of explicitly calculating 
conditional Hamiltonians depending on the block spin variables $ on coarser layers 
Afc. This formulation is possible if the conditional Hamiltonians are of the same 
type or similar to the fundamental Hamiltonian. Then, the conditional probabilities 
used for the updating on the k-th layer can be computed without always going 
back to the finest level Aq. Therefore, an explicit multigrid implementation reduces 
the computational work on the coarser layers (see the work estimates below). At 
least in free field theory, an explicit multigrid implementation is possible using 9- 
point prolongation kernels in two dimensions and generalizations thereof in higher 
dimensions |[T^, [T^. Generally, an explicit multigrid implementation for interacting 
models is only feasable in special cases with piecewise constant kernels. 

An algorithm formulated in the explicit multigrid style can always be translated 
to the unigrid language (that is how we are going to use the unigrid formulation). 
The reverse is not true, since not all nonlocal changes of the fundamental field 
configuration can be interpreted as updates of a single block spin variable of an 
explicit multigrid. As an example, one can use overlapping blocks in the unigrid 
style by translating the fields by a randomly chosen distance 

If we formulate our kinematical analysis in the unigrid language we nevertheless 
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can include all algorithms formulated in the explicit multigrid style. 

The sequence of sweeps through the different layers Ak of the multigrid is or- 
ganized in a periodic scheme called cycle [|T^. The simplest scheme is the V-cycle: 
The sequence of layers visited in turn is Aq, Ai, . . .Ak,Ak-i ■ ■ • Ai. More general 
cycles are characterized by the cycle control parameter 7. The rule is that from an 
intermediate layer A^ one proceeds 7 times to the next coarser layer A^+i before 
going back to the next finer layer Ak-i. A cycle control parameter 7 > 1 samples 
coarser layers more often than finer layers. With 7 = 1 we obtain the V-cycle. 7 = 2 
yields the W-cycle that is frequently used with piecewise constant kernels. 

The computational work estimates for the different cycles dbTQ clS follows 
The work for an explicit multigrid cycle is ~ L*^ if 7 < l'^, where L denotes the 
lattice size. The work for a unigrid cycle is ~ L'^logL if 7 = 1, and ~ L'^+'^^^n if 
7 > 1. Here, / denotes the blocking factor used in the iterative definition of the 
block lattices. 

If one wants the computational work in the unigrid style to not exceed (up to a 
logarithm) an amount proportional to the volume L'^ of the lattice, one has to use 
a V-cycle. Simulations with 7 > 1 (e.g. a W-cycle) can only be performed in the 
explicit multigrid style. 



3 An approximation formula for 

In this section we shall derive an approximate formula for the quantity Q{s) defined 
in (P). We can write Q{s) as 

n{s) = Jdu min(l, e"") J ^ e-'^" (e'^'^^) . (12) 

Let us assume that the probability distribution of A7i is approximately Gaussian. 
We parameterize this distribution as follows: 

dprob(AH) oc dAH exp(--^(A7^ - hiY) , (13) 

with hi = (AH) and = UiAH"^) - (AH)^). Then 

{ex-p(ipAH)) ~ ex-p{ihip — /i2P^) • (14) 

The integrations in eq. (|T2D can be performed exactly since there are only Gaussian 
integrals involved. The result is 

nis) ^ i (erfc(^) +exp(/.2-/.i)erfc(^^^)) , (15) 

with erfc(x) = 2/y^ I^dt exp(— t^). We shall now exploit the translational invari- 
ance of the measure Vcj) = V{(f) + sip) to show that the difference of hi and /i2 is of 
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order s^. The starting point is the observation that (exp(— A7i)) = 1. This imphes 
that 



^lu(exp(-AH)) = i: ((-AHD. =0. (16) 

s=0 m=l s=0 

((.))c denotes the connected (truncated) expectation value. Note that there are no 
contributions in the sum for m > n. This follows from the fact that A7i is of order 
s. Consequently, (A7i)™ = 0{s"^), and all contributions with m > n vanish in the 
limit s — > 0. For n=2 we obtain the relation 



UAH-). 



s=0 



ho 



0. 



[11] 



s=0 



If we assume that hi and /12 are even in s (which is the case if Ti. is even in 0) , then 
eq. (|T7|) says that the difference of hi and h2 is of order s^. 

We shall later demonstrate that the approximation hi ^ /i2 is in practice very 
good, even for small block size. In this case the acceptance rate prediction simplifies 
further, 

fi(s) ^ eiic{\^i) . (18) 
(For an analog result in the context of hybrid Monte Carlo see ||I9[ . ) 

For free massless field theory with Hamiltonian 7i(0) = ^(0, — A0), we get hi = 
^2 = with a = — Alp), and our approximation formula becomes exact 



Q{s) = erfc 



Eq. (|T9D can be checked directly by using {exp{ipAH)) 
(|r^. This relation is exact in free field theory. 



(19) 

exp{ihip — /i2P^) in eq. 



4 Coarse-to-fine interpolation 

In this section we shall discuss several choices of the coarse-to-fine interpolation 
kernels. In order to have a "fair" comparison, all kernels ip will be normalized 
according to eq. (H). 

In free massless field theory, the quantity a = [tp, —Aip) characterizes the de- 
crease of the acceptance rate VL{s) eq. (11^ with increasing shift s. Therefore it is 
natural to minimize a in order to maximize VL{s) for fixed s. 

The optimal kernel ip'^^"-^* from the point of view of acceptance rates can be 
defined as follows: minimize the quadratic form 

a = i^p,-Ai;) (20) 

under the constraints that the average of ip over the "central block" x'^ is given by 
^g-d)/2^ and its average over blocks x' 7^ x'^ vanishes: 

Lb' E = I.g''^/'<5.',.^ for all x' e Ak . (21) 
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This variational problem can be solved with the help of Fourier methods. The result 
is 

2+d 

rr^ = , (22) 

where Ax,x'^ denotes the Gaw§dzki-Kupiainen kernel (see, e.g. 0). The use of this 
kernel leads to a complete decoupling of the different layers of the multigrid. This 
way of interpolating from a coarser block lattice to the fine lattice Aq is well known 
in rigorous renormalization group theory [^. It is interesting that considerations 



about optimizing acceptance rates in a stochastic multigrid procedure lead to the 
same choice of the interpolation kernel. 

Because ip'^^"-'^^ is nonvanishing on the whole lattice, it is not convenient for 
numerical purposes. For an attempt to change the block spin $2,/, on block x'^ one 
has contributions to the change of the Hamiltonian from all lattice points. Therefore 
the computational work for a single update is proportional to the volume. 

We define a "truncated kernel" 7^*™"^ by restricting the support of tp on the 
block x'^ and its nearest neighbor blocks y'^ 

^trunc ^ Q jf X ^ x'^ or X ^ y'^, where y'^ n.n. x'^ . (23) 

In other words, the Laplacian in eq. ( ^0|) is replaced by a Laplacian with Dirichlet 
boundary conditions on the boundary of the support of ip. We again minimize 
a = [ip, —Antp) under the 2d + 1 constraints that the average of ip over the blocks 
Xg and its nearest neighbor blocks is given. This minimization can be performed 
numerically by a relaxation procedure. In order to maintain the normalization 
condition, one always updates simultaneously two spins residing in the same block, 
keeping their sum fixed. The '?/;*'^""'^-kernels were used in a multigrid simulation of 
the 0^ model in four dimensions ^1 . 



From a practical point of view, it is convenient to use kernels that have support 
on a single block x'^, i.e. 

= Oifx^x;. (24) 

We define a kernel ip"^"^^ with this property by minimizing a = {ip, —Ad^^iP) under 
the constraint that the average of ip over the block x'^ is given. The Laplacian with 
Dirichlet boundary conditions on the boundary of x'^ is defined as follows: 



-2d(Px+ E 



y n.n.x 



for X & x'g . (25) 



^mm ^g^j^ calculated using an orthonormal set of eigenfunctions of A/),^'^- 

We shall now discuss other kernels with support on the block that are frequently 
used in the literature. 



piecewise constant interpolation: 



i.const J Lb^"^ '^^^'^ for X G 



^const ^ J - - (26) 

10 for s ^ Xg . 
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Table 1: Results for a = (\E', — A\E') in 2 dimensions, 512^ lattice 



kernel 


Lr=2 


Lr=4 


Lr=8 


Lr=16 


Lb=32 


Lr=64 


Lb=128 


Lr=256 


exact 


C QCiCi 

D.oyy 


o.yuz 


y. / UO 


y.y4i 


iU.UU 


1 n no 


iU.io 


io.ii 


trunc 


7 nnn 
( .UUU 


o /I nc; 
y.4U0 


iU. ( o 


ii.oo 


ii.oy 


ii.o4 


1 1 oo 

i i .yz 




min 


8.000 


13.24 


18.48 


22.58 


25.23 


26.76 


27.59 


28.02 


sine 


8.000 


13.62 


19.34 


23.78 


26.62 


28.25 


29.13 


29.58 


linear 


8.000 


15.80 


24.58 


31.84 


36.68 


39.51 


41.05 


41.84 


const 


8.000 


16.00 


32.00 


64.00 


128.0 


256.0 


512.0 


1024 



This kernel has the advantage that for many models the conditional Hamiltonians 
used for updating on coarse grids are of the same type or similar to the fundamental 
Hamiltonian. This means that the conditional probabilities used for the updating 
on the k-th layer can be computed without always going back to the finest level Aq. 
Therefore, an explicit multigrid implementation with a W-cycle can be used. 



piecewise linear interpolation: 
We consider the block 



x'^ = {x\x^ G {1, 2, 3, ... , Lb}, fi = I, . . . ,d} . 



(27) 



The kernels for other blocks are simply obtained by translation. For Lb even, 7^'*"^'*'' 
is given by 



^linear = I + ^ 

11=1 

is a normalization constant 



Lb + 1 



for X G x' 



{2i 



ground state projection kernels: 

^sme ^YiQ eigenfunction corresponding to the lowest eigenvalue of the negative 
Laplacian with Dirichlet boundary conditions — A^^^./^ : 



TT 



sm 



'Lb + 1 







for X G x' 



for X ^ x'g 



(29) 



Again, A/" denotes a normalization constant. Note that this kernel is different from 
^m«n_ js^ generalization of this kernel was introduced for scalar fields in the back- 



ground of nonabelian gauge fields in ref. ||2^ . 
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Table 2: Results for a = — A\l/) in 4 dimensions, 64'^ lattice 



kernel 


Lr=2 




Lr=8 


Lr=16 


Lb=32 


exact 


T /I /I Q 


on Qo 


OQ /I Q 


z4. / 1 


on f^o 
oU.DZ 


trunc 


14.0 ( 


/l.ul 


ZD. 04 


OO OR 

zy .ZD 




min 


16.00 


27.72 


41.56 


54.18 


63.33 


sine 


16.00 


30.37 


48.46 


64.85 


76.44 


linear 


16.00 


39.02 


70.78 


101.0 


122.9 


const 


16.00 


32.00 


64.00 


128.0 


256.0 



The results for the quantities a = {ip, —^ip) fo^' different kernels in two dimen- 
sions are presented in table |l|. We used a 512^ lattice (t^^^'^^* depends on the lattice 
size). The different kernels are ordered according to increasing value of a. 

The values of a^^'^'^* and a*™"^ are close together. This shows that the truncation 
of the support of ip to the block and its nearest neighbor blocks is a good approxi- 
mation to (in the sense of acceptance rates). The value of a^^'*^* for Lb = 256 
is remarkably higher than on smaller blocks. This is a finite size effect because the 
block lattice consists only of 2^ points. Since the nearest neighbors overlap on a 2^ 
lattice, no result for q/*'^""'= is quoted for Lb = 256. The values of a for the smooth 
kernels with support on the block ^^"-^ ^sme ^Unear same magnitude. 
We can see that ip'^^^'^ is almost as good as the optimal ip = ip"^^^. 

The results for different kernels in four dimensions are presented in table Here 
we used a 64*^ lattice. In principle, the a's behave as in two dimensions. The values 
of tt'*"^'^'' for small blocks are higher than a™"**. The pyramids of the piecewise 
linear kernels have a lot of edges in four dimensions which lead to high costs in the 
kinetic energy. 

The L^-dependence of the a's in d dimensions is 



a = 2dLB for piecewise constant kernels . 
a — > const for smooth kernels . 

Lb»1 

As an example, the expression for a''*"'^ in d dimensions is 



(30) 



Ll+''{LB + lY2''+^dsm"'+^ 
For large block sizes we find 



TT 



2(Lb + 1) 



a'^ne — ^ d—— = const . (32) 

Lb»1 23rf ^ ^ 

From table |l] we observe that in two dimensions a becomes almost independent of 
Lb for the smooth kernels if the block size is larger than 16. In four dimensions 
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(table 1^), we find a{LB) ~ const only for a^^'"'^*. The other a's for the smooth 
kernels have not become independent oi Lb for the block sizes studied. 



5 Acceptance rates in free field theory 



IS 



Recall that we have VL{s) = erfc \s\j in massless free field theory. Q{s 

only a function of the product as"^. In order to keep Q{s) fixed when the block size 
Lb increases we have to rescale the changes s like ^(Lb)^^/^. As a consequence, 
to maintain a constant acceptance rate in massless free field theory, s has to be 
scaled down like Lb for piecewise constant kernels, whereas for smooth kernels 
the acceptance rates for large do not depend on the block size. 

Note that this behavior of the acceptance rates for large is not yet reached in 
four dimensions for the block sizes studied (except for ■?^^^'*^*) . See also the discussion 
of the Metropolis step size below. At least for free field theory, the disadvantage of 
the piecewise constant kernels can be compensated for by using a W-cycle instead of 
a V-cycle. Smooth kernels can be used only in V-cycle algorithms. An exception are 
9-point prolongation kernels in two dimensions and generalizations thereof in higher 
dimensions. They can also be used with a W-cycle, at least in free field theory (cf. 
section |^) . 

We now illustrate what this rescaling of s means for the Metropolis step size e 
in an actual multigrid Monte Carlo simulation. Look at the integrated acceptance 
probability defined in eq. (pT]). If we insert the exact result (p!9| ) for massless free 
field theory we get 

(e) = erfc ( ^f^e) + \l - e-f^'l . (33) 



P.. 



Pace is only a function of the product ae^. In order to keep Pace fixed (to, e.g. 50 
percent) we have to rescale e{LB) like a(LB)~^/^, exactly in the same way as we had 
to rescale s to keep Q{s) fixed. This L^-dependence is plotted in figure 1 for two 
dimensions and in figure 2 for four dimensions. 



We now discuss massive free field theory with Hamiltonian 7i(0) 



m^]0). We find hi 



1 



Or 



2 



, s , with am given by 

= {ip, [-A + m^]i!) =a + m^ ^ ipl 

a;GAo 



\ [-A + 



(34) 



Therefore the exact result is Vl[s) = erfc (^^Jam/S \s\j. A term J^xi'l scales ~ in 
arbitrary dimensions. For piecewise constant kernels 



const 



For '?/''**"'^-kernels we find 
, 2 



a;GAo 



sme 

X 



2<i^2+c( 



xGAo 



TT 



2{L 



B 



sm 



-2d 



TT 



Lf 



(35) 



(36) 
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and for large block sizes 



9 -7r2d 



If the block size L_b is smaller than the correlation length C, = 1/m, hi is still 
dominated by the kinetic term (■?/', —Aip), and the discussion is the same as in the 
massless case. 

As soon as the block size Lb becomes larger than ^, hi is dominated by the mass 
term s^m^ Y,x i^l. ~ ■^^L^, and s has to be rescaled like s ~ L~j^ in order to maintain 
constant acceptance rates. Of course this is a dramatic decrease for large block sizes 
compared to s ~ const (using smooth kernels) in the massless case. Block spins on 
large blocks are essentially "frozen" . But this is not dangerous for the performance 
of the algorithm in massive free field theory: The effective probability distribution 
for the block spins $ is given by exp(— 7Yofr('^')), where ?ieff(^) denotes the effective 



Hamiltonian in the sense of the block spin renormalization group . The physical 



fluctuations of the block spins are dictated by an effective mass term 

mis E ^l' with mil, ~ m'Ll . (38) 

Thus, the algorithmic fluctuations (described by the mass term rn? ''Pi ~ m?L\) 
and the physical fluctuations (described by the effective mass ~ m'^L'^) behave 
similar, and the multigrid algorithm is able to create fluctuations just of the size 
that is needed by the physics of the model. Moreover there is no need to do updates 
at length scales larger than ^ in order to beat CSD. 

In this sense, the discussed algorithmic mass term J2x is well behaved for 
free field theory, since it decreases with the physical mass in the vicinity of the 
critical point. As we shall see in section for interacting models close to criticality, 
a different scenario is possible. There, it can happen that an algorithmic "mass 
term" ~ Y^x'^Pl persists, whereas the renormalized mass vanishes. If this happens, 
the multigrid algorithm is not able to produce the large critical fluctuations required 
by the physics, and we can not expect that CSD will be eliminated. 

The L^-dependence of a term J^xi^t ^^^^ ^l^o be needed in the study of the 0^ 
theory in section]^ below. In d dimensions such a term scales ~ L^"^: For piecewise 
constant kernels 

E {rr^'Y = , (39) 
whereas using ?/'**"''^-kernels we find 



j2 (^--*)' = 6'^L^+2"(LB + l)"sin« 

xGAq 



IT 



.2(Lb + 1) 

In the limit of large block sizes this term behaves like 



xgAo 
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In order to summarize the different large-L^-behavior of local operators in tlie 
kernel ip discussed here, let us introduce the degree of relevance in the sense of the 
perturbative renormalization group: The (superficial) degree r of relevance of a local 
operator in ip which is a polynomial of m scalar fields with n derivatives is defined 
hj r = d + m{2 — d)/2 — n. This definition is valid for smooth kernels. For large 
Lb, an operator with degree of relevance r behaves like L^. An operator is called 
relevant if r > 0. As we have seen in the examples above, a mass term has r = 2, 
and a ip^-term has r = A — d. A kinetic term a = {ip, —Aip) has r = for smooth 
kernels. 

The only difference for piecewise constant kernels is that a kinetic term behaves 
like a = {ip, —Aip) oc Lb- 



6 Acceptance rates for interacting models 

In this section, we shall apply formula (|r^) in the discussion of multigrid procedures 
for three different spin models in two dimensions: the Sine Gordon model, the XY 
model, and the single-component 0^ theory. The scale dependence of acceptance 
rates for interacting models will be compared with the behavior in free field theory, 
where CSD is known to be eliminated by a multigrid algorithm. 

6.1 2-dimensional Sine Gordon model 

The 2-dimensional Sine Gordon model is defined by the Hamiltonian 

^(0) = ^(0, -A0) - CE . (42) 

The model undergoes a Kosterlitz-Thouless phase transition at jSc- In the limit 
of vanishing fugacity (, the location of the critical (3 is exactly known: (3c Svr for 
C — > 0. For (3 > (3ci the model is in the massless phase, and the flow of the effective 
Hamiltonian (in the sense of the block spin renormalization group) converges to that 
of a massless free field theory: the long distance behavior of the theory is that of a 
Gaussian model. Since multigrid algorithms have proven to be efficient in generating 
long wavelength Gaussian modes, one might naively conclude that multigrid should 
be the right method to fight CSD in the simulation of the Sine Gordon model in the 
massless phase. But this is not so. For hi we find the expression 

/^i = ^^' + CCE[l-cos(#.)], (43) 

with C = (cos 0a;). Recall that hi is the quantity that determines the acceptance 
rates Q{s): 

n{s) ^ eikil/hi) . (44) 

The essential point is that the second term in ( ^3|) is proportional to the block volume 
for piecewise constant and for smooth kernels (cf. the discussion in section 



12 



This can be checked for small s by expanding in s. One therefore has to face 
a dramatic decrease of acceptance when the blocks become large, even for small 
fugacity (. A constant acceptance rate is achieved only when the proposed steps are 
scaled down like L~^^. It is therefore unlikely that any multigrid algorithm - based 
on nonlocal updates of the type discussed in this paper - will be successful for this 
model. 

We demonstrate the validity of formula (^) (using a Monte Carlo estimate for 
C) by comparing with Monte Carlo results at /3 = 39.478 and ( = 1- This point is 
in the massless phase, where the correlation length ^ is or the order of the lattice 
size L. In figure 3 we show both the numerical and analytical results for Q{s) for 
Lb = 4, 8, 16, 32 on lattices of size 16^, 32^, 64^, 128^, respectively. 

We tested the precision of our approximation formula for piecewise constant 
kernels only. However, we have no doubts that the quality of the approximation is 
also very good for other shape functions ip. 



6.2 2-dimensional XY model 

We now discuss the 2-dimensional XY model, defined by the partition function 

z= fUde, exp(/5 cos(e, - e,)) . (45) 

X <x,y> 

The sum is over all unordered pairs of nearest neighbors in the lattice. As the Sine 
Gordon model, the XY model has a massless (spin wave) phase for (3 > (3c-, and a 
massive phase for j3 < (3c- The best available estimate for the critical coupling is 
/3e = 1.1197(5) m. 

Nonlocal updates are defined by 

Qx^Qx + s^x , (46) 
with ip obeying again the normalization condition (§). To define a (linear) block 



spin, we rewrite the partition function (|45|) in terms of 2-component unit vector spin 
variables Sx'- 

Z= fU {£s. 5isl - 1)) expiP ■ ^y) (47) 

X <x,y> 

The block spins S^' are then defined as block averages of the unit vectors s^- 

Note that the proposal ( |^ ) changes the block spin by an amount ~ s only when 
the spins inside the block are sufficiently aligned. This will be the case in the spin 
wave phase for large enough (3- For smaller (3, the correct (or "fair") normalization 
of the kernels ip is a subtle point. We believe, however, that our argument is not 
affected by this in a qualitative way. 

The relevant quantity hi is given by 

/zi = ^ [1 - cos(s(^, - ^y))] , (48) 

<x,y> 
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with E = (cos(G2: — Oy)), x and y nearest neighbors. For piecewise constant kernels, 
hi is proportional to Lb- For smooth kernels hi will become independent of Lb for 
large enough blocks. For small s, 

hi ^ \s^(3E (^x - i^y? = Wf^Ea . (49) 

<x,y> 

As above, a = (■?/', —Aip). This quantity becomes nearly independent of already 
for Lb larger than 16 (cf. section 

From the point of view of acceptance rates the XY model therefore behaves like 
massless free field theory. A dynamical critical exponent z consistent with zero was 
observed in the massless phase ||2^. The failure of multigrid Monte Carlo in the 



massive phase {z ~ 1.4 for piecewise constant kernels ||25|) is an example for the 
fact that good acceptance rates are not sufficient to overcome CSD. 

We again checked the accuracy of formula (|TBp by comparing with Monte Carlo 
results at /5 = 1.2 (which is in the massless phase, where the correlation length ^ 
is of the order of the lattice size L). The only numerical input for the analytical 
formula was the link expectation value E. The results are displayed in figure 4. 

One can do a similar discussion for the 0{N) nonlinear cx-model with N > 2, 
leading to the same prediction for the scale dependence of the acceptance rates. This 
behavior was already observed in multigrid Monte Carlo simulations of the 0(3) 
nonlinear a-model in two dimensions with smooth and piecewise constant kernels 
12611. 



6.3 2-dimensional 0^ theory 

Let us now turn to single- component ci-dimensional 0^ theory, defined by the Hamil- 
tonian 

X ^- X 

For hi one finds 

/^i = + + ^p] E^^} + ^ E v^^ (51) 

where P = {(pl). We have used that expectation values of operators which are odd 
in (f) vanish on finite lattices. Recall that J^xi^x increases with L^, independent of 
d, whereas Y.x i't scales like L'^'^, for smooth and for piecewise constant kernels (cf. 
the discussion of the different kernels in section We conclude that also in this 
model we have to face rapidly decreasing acceptance rates when the blocks become 
large. As in the case of the Sine Gordon model, s has to be rescaled like in 
order to maintain constant acceptance rates. 

Therefore there is little hope that multigrid algorithms of the type considered 
here can overcome CSD in the 1-component 0"^ model. In numerical simulations of 2- 
dimensional 0^ theory, a dynamical critical behavior is found that is consistent with 
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2 for piecewise constant and for smooth kernels H O, p^. In four dimensions, 



there is no definite estimate for z [21 



Figure 5 shows a comparison of Monte Carlo results for 2-dimensional 0^ 
ory with the theoretical prediction based on the numerical evaluation of P. 
simulations were done in the symmetric phase at = —0.56 and Ao = 2.4. 
correlation length at this point is ^ ~ 15 [B^ . 



the- 
The 
The 



6.4 Summary of section 6 

Our approximation formula has proven to be quite precise. The results for three 
different models are consistent with the following rule: 

Sufficiently high acceptance rates for a complete elimination of CSD can only be 
expected if hi = {H^cj) + sip) — H^cp)) contains no relevant operator in ip. 

As we have seen above, the typical candidate for a relevant operator in ijj is 
always an "algorithmic mass term" of the type ~ J2x with degree of relevance 
r = 2. 

This rule is formulated for smooth kernels. For piecewise constant kernels, it has 
to be modified. There, the kinetic term a oc is relevant as well. At least in free 
field theory this disadvantage can be compensated for by using a W-cycle. Apart 
from this modification the rule carries over to the case of piecewise constant kernels. 



7 A multigrid procedure for lattice gauge fields 

In this section we propose a multigrid procedure for pure lattice gauge theory and 
study the behavior of acceptance rates with increasing block size Lb- 

We consider partition functions 

Z= flldU,^,exp{-nm)- (52) 

X,/i 

The link variables U^:^^ take values in the gauge group U{1) or SU{N), and dU 
denotes the corresponding invariant Haar measure. The standard Wilson action 
Hi^U) is given by 

= /^Efl- jr^eTiUr] . (53) 

V 

The sum in ( ^3|) is over all plaquettes in the lattice. The U-p are path ordered 
products around plaquettes V, 

Ur = U,,,U.+^,MU,^,Ul, . (54) 

U* denotes the hermitean conjugate (= inverse) of U. 
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7.1 The abelian case 



We now consider the case of gauge group U{1). The hnk variables are parameterized 
with an angle — vr < O^^fj, < vr through 

[4,^ = exp(z6'a:,^) . (55) 

Nonlocal updates can be defined as follows: One chooses a hypercubic block x'^ of 
size and a direction r with 1 < t < d. During the update, r will be kept fixed. 
All the link variables Ux^r attached to sites x inside the block x'^ are proposed to be 
changed simultaneously: 

f4,r — > exp{isilJx)Ux,r ■ (56) 
Again, ip denotes an interpolation kernel as introduced in section ^. This updating 



scheme was introduced and studied in two dimensions in ref. 131 and also in four 



dimensions in ref. |]T3, using piecewise constant kernels. Of course, one can use all 



versions of smooth kernels as well, with their support not necessarily restricted to 
the block x'^. 

Let us now study the acceptance rates for these update proposals with the help 
of formula (pisf ). We consider general kernels ip. For hi = (A7i) we find 

^1 = /^^ E E [1 - cos(s(V^,+^ - ij,))] , (57) 

with P = (TtU-p). The sum does not include contributions from links which point in 
the fixed r-direction. If we denote the "slice" of lattice sites with r-component t as 
A[ = {a; G Ao I Xt = t}, we see that hi is a sum of independent contributions from 
slices orthogonal to the r-direction. Therefore, no smoothness of kernels is needed 
in the r-direction, and from now on we choose V'x to be constant in that direction. 
Let us assume that the support of ip in r-direction is contained in the block x'^. 
Then we find for small s 

hi ^ \s^(3PLb E E(^-+A - V-x)' = Is'fiPad-i . (58) 
with ad-i = , —Alp'). Here, ijj' denotes the kernel ip restricted to the d — 1 

1/2 

dimensional sublattice A[, multiplied with a factor in order to be properly 
normalized as a (c? — 1) dimensional kernel. 

From the kinematical point of view, the behavior of acceptance rates in the U{1) 
lattice gauge theory in d dimensions is the same as in massless free field theory. One 
might therefore expect that it is possible for a multigrid algorithm to overcome CSD 
in this model. Indeed, in numerical simulations in two dimensions using piecewise 
constant kernels, the dynamical critical exponent was found to be 2; ~ 0.1 p!3| . 



However, it was also observed that the multigrid algorithm is not able to move 
efficiently between different topological sectors. The above quoted exponent should 
therefore be interpreted with some care. For the results of a study of a multigrid 
algorithm for 4-dimensional ?7(l)-theory see . 
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Let us conclude the discussion of the abehan case with the remark that with 
no loss of generality one could consider blocks x'^ that consist only of one layer in 
r-direction, i.e., effectively (ci — l)-dimensional blocks. This is a consequence of the 
fact that the updates of the two variables Ux,t and Vy^.^ are statistically independent 
if Xt 7^ yr- This property carrys over to the nonabelian case. 

7.2 The nonabelian case: gauge group SU{2) 
7.2.1 Covariant nonlocal update proposal 

We shall now discuss a generalization of the above described procedure to the non- 
abelian case. To be concrete, we study 4-dimensional SU (2) lattice gauge theory. 

Let us start with an attempt of a straightforward generalization of the procedure 
described for the abelian theory. This would amount to propose updates 

f/,,.-.f/;, = i?, (59) 

where the "rotation" matrices are in SU(2). We parametrize them as 

Rx{n,s) = cos^sipx/'^) + isin^sipx/^) n-cr , (60) 

where n denotes a three-dimensional real unit vector, and the ai are Pauli matrices. 
ijj will have support on 3-dimensional blocks x' of size L%, and the blocks consist of 
points lying entirely in A[. 

We use the fact that updates of link variables in different slices are statistically 
independent (as discussed at the end of subsection |7l|). One possible updating 
scheme is to perform the proposed updates on different slices in sequence. Another 
possible updating scheme consists of building hypercubic four-dimensional blocks out 
of "staples" of Lb three-dimensional blocks of size and to perform the updates 
on this hypercubic block simultaneously. Because of the independence of different 
slices, the analysis of acceptance rates is the same for both cases. For simplicity we 
study three-dimensional blocks here. 

The energy change associated with the update proposal (^) is 

= -f E T^^iUv - t^p) = -f E E Tr{iR:Ux,,Rx+i. - Ux,,)H:J , (61) 

with H*^^ = Ux+fi^TU*^^^^U* ^. The relevant quantity for the acceptance rates is 
hi = (AH) . For piecewise constant kernels ip one gets 

hi = ^A' {Lb - l)Ll sm\sL~s'^'/2) + 3/5PL|[1 - cos(sL^'/V2)] , (62) 

with 

A' = -{TTiin-aUx,,n-a-Ux,,)Hl^)), 

P = {Tt{Ux,,H:J) = (TrUr) . (63) 
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The first contribution to A' is the expectation value of a quantity that is not 
gauge invariant. Determining its gauge invariant projection, we can show that this 
contribution vanishes: 

{Tiia-nU,,^a-nHl^)) = J dV{TT{n-aVU,,,n-aHl^V*)) = 0, (64) 
because for SU{2) 

J dVTriAVBV") = ^TtATtB , (65) 

and the Pauh matrices are traceless. Therefore we get A' = P. 

To the first term in eq. ( |62D all links contribute that are entirely inside the block, 
whereas the second term sums the contributions of all links that have one site in 
common with the surface of the block. For small s, the first term behaves like s^L^. 
This is exactly the behavior of a mass term that, as we have learned in the previous 
sections, is toxic for the multigrid algorithm. 

The main difference to the abelian case is that in addition to the costs from the 
surface of the block we have a contribution from the interior of the block. Unfortu- 
nately, this contribution grows quadratic with the block dimension Lb- 

This does not come as a surprise. Due to the gauge invariance of the model, 
the L'^'s do not have a gauge invariant meaning. Therefore the rotations Rx that 
are smooth over the block for a given gauge can be arbitrarily rough after a gauge 
transformation. It is therefore clear that the rotation matrices have to be chosen in 
a gauge covariant way. 

We generalize the update proposal (|59D as follows: 

Ux,T ^^,r = 9x,Rx9x Ux,T , (66) 

with Qx G SU{2). Note that in the abelian case we obtain nothing new, because Qx 
and Rx commute. In the nonabelian case, we find for piecewise constant 

/^i = f E sm'isL-^'^'/2) + SPPLIII - cos(.L^^/V2)] , (67) 

{x,x+fi)ex'^ 

with 

Ax,, = -{Tiiin.aUl^n.a - f/f,,)iff;)) . (68) 

Here we have introduced the notation U^^^ = gxUx,f_igx+fi- is defined analogously. 

is the gauge field obtained by applying a gauge transformation g to U . We are 
free to choose g. To obtain a valid Monte Carlo algorithm, we require that the 
g^s should not depend on the link variables to be updated, i.e. those living on the 
links (a;,a; + f). On the other hand we want to minimize hi. Let us inspect the 
quantity Ax,, defined in eq. (|68|) that leads to the unwanted mass term behavior of 
hi. Consider the extreme case of — >■ cx). Then the allowed configurations are pure 
gauges, i.e. configurations that are gauge equivalent to Ux,, = 1 for all x, /x. If we 
choose g as the transformation that transforms all links to unity, it is obvious that 
Ax,, is zero. In particular, to obtain this result, it is sufficient to gauge all links 
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inside the block that do not point in the r-direction to unity. This consideration 
leads to following proposal: Choose g as the gauge transformation that maximizes 
the functional 

Gc,.',{U,g)= E Tr{g,U,,,g*^^) . (69) 

We call this gauge "block Coulomb gauge". This gauge will bring the links inside 
the block as close to unity as possible thus leading to a kind of minimization of A^^^ 
(corresponding to a minimization of the mass term). Note however, that we do not 
intend to actually perform the gauge transformation. We use the concept of gauging 
only to define covariant rotations g^RxQ^- Covariancc here means that the relevant 
quantity Ti{{n-aU^^^n-(j — U^^^)H^*^ is now gauge invariant. To see this, assume 
that we pass from U to by applying the gauge transformation h. The Coulomb 
gauge condition will then lead to a new g' — gh*. Now note that [U^y^* — . The 
same argument applies to H . 

Let us summarize the steps of the nonlocal updating scheme for SU{2): 

1. Choose a block x'^ of size L\ that is contained in the shce A[. All hnk vari- 
ables Ux,T pointing from sites x inside the block in r-direction will be moved 
simultaneously. 

2. Find the gauge transformation g defined by the block Coulomb gauge condition 

Gc,x',{U,g)= ^(^xC/x,M^:+/i) = maximal. (70) 

{x,x+ij,)&'„ 

3. Propose new link variables C/^^ by 

U^,r^K,r = glRxgxU^,r, (71) 

with 

Rx{n, s) = cos{sipx/'^) + i smi^sTpx/'^) n-a . (72) 

s is a uniformly distributed random number from the interval [—e, e], and n is 
a vector selected randomly from the three dimensional unit sphere. 

4. Calculate the associated change of the Hamiltonian ATC and accept the pro- 
posed link variables with probability min[l, exp(— A?i)]. 

The detailed balance condition is fulfilled by this updating scheme: For the 
naive version with = 1 it is straightforward to show that the detailed balance 
condition holds, since the rotation matrices are chosen according to a probability 
distribution which is symmetric around unity. 

If we now take g according to some gauge condition, we have to be careful that 
we get the same g before and after the move Ux,t U!^,t- Otherwise this move would 
not be reversible. In other words: The gauge condition yielding g must not depend 
on Ux^T- This is indeed the case, since only link variables Ux,^ with [i ^ t enter in 
the block Coulomb gauge functional. Note that we do not have to fix the gauge 
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perfectly. If we always use the same procedure in finding g (e.g. a given number 
of relaxation sweeps starting from g = 1), we will always get the same g and the 
nonlocal update is reversible. 

As usual we now choose different (possibly overlapping) blocks x', different block 
sizes Lb, different slices A[ and different orientations r of the slices in turn. 

7.2.2 Acceptance analysis for nonlocal S'f/(2)-updates 

First numerical studies revealed that there is no substantial difference in the accep- 
tance rates when instead of using the block Coulomb gauge condition one uses the 
Coulomb gauge condition for the whole slice A[: 

Gc{U,g)= J2 Tr(5(^?7^,^5(*^^) = maximal. (73) 

{x,x+fi)£Al 

From a practical point of view this gauge condition is very convenient, because the 
relaxation algorithm to determine the gx can then be vectorized in a straightforward 
way. 

If we use the gauge condition (0), the quantity A^^^ becomes translation invari- 
ant and also independent of /i (where we still keep jj, ^ t). We get 

h, = {Lb - l)Ll sm\sL-^'/'/2) + 3/3PL|[l - cosisL-^'^' /2)] , (74) 

with 

A = -(TA{n-aUl^n.a - Ul^)Hi:^)) (75) 

Following the discussion after eq. (^), we identify the square root of A with a 
"disorder mass" m^i, 

mD = ^- (76) 

rriD has the dimension of a mass. It vanishes in the limit /? ^ oo, just like 
physical masses in the theory. Because of the disorder inside the block, m^) is 
nonzero for finite (3. This would not be a problem if for large correlation length tud 
scaled like a physical mass (cf. the discussion for massive free field theory in section 

I)- 

7.2.3 Monte Carlo study of m^j 

We computed m£, hy Monte Carlo simulations for several values of (3. To maxi- 
mize Gc we used a standard Gauss-Seidel relaxation algorithm vectorized over a 
checkerboard structure. The relaxation procedure consists in going through the lat- 
tice and minimizing the gauge functional ([73|) locally. For production runs it would 
be advantageous to use an accelerated gauge fixing algorithm such as overrelaxation 
or multigrid In the Monte Carlo studies reported in this section, we always 

used 50 GauB-Seidel sweeps to determine g. Note that by this procedure, Gc is not 
entirely maximized, especially on very large lattices where the relaxation algorithm 
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Table 3: Monte Carlo results for rriD and P 



lattice size 




rriD 


P 


statistics 


8^ 


2.4 


0.507(2) 


1.5131(6) 


10,000 


12^ 


2.4 


0.4957(4) 


1.5121(3) 


10,000 


16^ 


2.4 


0.4955(2) 


1.5119(1) 


10,000 


8^ 


2.6 


0.497(4) 


1.7429(3) 


30,000 


12^ 


2.6 


0.465(2) 


1.7424(2) 


20,000 


16^ 


2.6 


0.4644(3) 


1.7421(1) 


10,000 


20^ 


2.6 


0.4650(2) 


1.7422(1) 


5,000 



Table 4: Comparison of m/j with physical masses 



lattice size 


P 


rriD 




mo+ 




m/)/mo+ 


16^ 


2.4 


0.4955(2) 


0.258(2) 


0.94(3) 


1.92 


0.53 


20'* 


2.6 


0.4650(2) 


0.125(4) 


0.52(3) 


3.72 


0.89 



suffers from CSD. However, for the detailed balance to be fulfilled, we only need that 
one uses always the same number of relaxation sweeps. Several tests revealed that 
increasing the number of relaxation sweeps beyond 50 did not affect the acceptance 
rates in a substantial way. In our implementation, 50 GauB-Seidel sweeps over all 
slices of a given direction r required the same CPU time (on a CRAY Y-MP) as 
four Creutz heatbath SU{2) update sweeps. 

We checked the validity of the acceptance formula (|I8|) using Monte Carlo esti- 
mates for m/5 and P. Figure 6 shows results for (3 = 2.6 on a 20^ lattice. The results 
perfectly justify the usage of the approximation formula. It is therefore sufficient 
to study the behavior of the quantities and P. Our Monte Carlo results are 
presented in table |[ The last column gives the statistics in sweeps (equilibration 
sweeps are not counted here). We used a mixture of four microcanonical overrelax- 
tion sweeps followed by a single Creutz heat bath sweep. Measurements (including 
the determination of g) were performed every 25 sweeps. 

In table ^ we display the ratios of the disorder mass m/j with two physical 
masses, the square root of the string tension k and the lowest glue ball mass mo+. 



The estimates for the physical masses are taken from ref. [^. The results show 
that the disorder mass is nearly independent of P in the range studied, whereas the 
physical masses decrease by roughly a factor of two. Thus, is not scaling like a 
physical mass for the couplings studied here. We conclude from this that for large 
blocks the term quadratic in will strongly suppress the acceptance rates even 
when the ratio of correlation length and block size Lb is kept constant. 
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However, one could hope that the value of the unwanted mass term is so small 
that it does no harm in practical calculations. Let us examine the effect of this 
mass term in more detail. Recall that hi is built up from two contributions. The 
first contribution is that related to the gauge field disorder inside the block and is 
quantitatively represented by the mass mr). The second contribution is associated 
with the block surface. The latter can of course be made smaller by using smooth 
kernels t/j instead of the piecewise constant kernels discussed so far. However,the 
disorder term cannot be expected to become smaller for smooth kernels (see below) . 
In figure 7 we plotted separately the two contributions to hi 

hi,A = ^A{LB-l)Llsm\sL-^'^'/2), 

hi,p = 3pPLl{l-cos{sLs'/y2)), (77) 

for P — 2.6 and block size = 8 on a 20^ lattice. The plot shows that already 
for this block size the disorder contribution is by no means small - it is comparable 
to the surface effect. It is therefore not clear that one could achieve any significant 
improvement by using smooth kernels. To investigate this in more detail, we derive 
an expression for hi, valid for smooth kernels as well: 

hi = ^s'A E + l:s\P - A)a, + O(s') . (78) 

Since Y^i^l ^ L\,we get essentially the same behavior for the disorder contribution 
as in the case of piecewise constant kernels. 

For smooth ■0**'*^ kernels we show separately in figure 8 the two contributions to 

hi 

hi,A = ^-s'A Y: , hi,p^A = ^s'{P - A)as , (79) 

for j3 = 2.6 and block size = 8 on a 20^ lattice. We observe that the surface 
effects are lowered by the smooth kernels, but the disorder contribution is even 
higher than for piecewise constant kernels. 

Piecewise constant kernels have the practical feature that once the change of the 
Hamiltonian has been calculated, one can perform multihit Metropolis updating or 
microcanonical overrelaxation. In a special case, even an explicit multigrid imple- 
mentation with a W-cycle is possible (sec below). For smooth kernels the change in 
the Hamiltonian would have to be calculated again and again. Also the advantage 
of smooth kernels are not that clear on small three or four dimensional blocks. For 
an actual simulation we would therefore prefer piecewise constant kernels. 

7.2.4 Mctximally abelian gauge 

Our proposal for the choice of g was motivated by the desire to minimize the quantity 
A. We now ask whether there is a better choice of g than the g determined by the 
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Coulomb gauge condition. For the sake of simplicity let us take n = (0,0, 1), i.e. 
n-a = (J3. Then A is given by 

A = -(Tr((a3 t/f,^ ^3 - f/,^,^)i/f;)) • (80) 

The choice of the Coulomb gauge condition aimed at bringing f/J'^ as close to unity 
as possible. Alternatively, one might require that f/|^ should be as close as possible 
to a S'f/(2)-matrix of the form + ia^a^. This will also lead to a small A. The 
corresponding gauge transformation g can be found by maximizing the functional 

GAiU,g)= 5: Tr(o-3f/,>3t^.';) , (81) 



leading to the maximally abelian gauge [^, here implemented only on a slice. We 



computed mo also using the g's resulting from this gauge condition and compared 
the results with the ones obtained by using the Coulomb gauge condition. We did 
not find a substantial difference. We prefer the Coulomb gauge condition because 
it does not depend on the direction n and thus saves computer time. 

7.2.5 Proposal for an implementation 

An explicit multigrid implementation is possible in a special case if we use piecewise 
constant kernels. This was pointed out in ref. [|I^ in a related context. 



The idea is to update only a fixed U{1) subgroup of SU{2) globally: We divide 
the fundamental lattice Aq in hypercubic blocks x' of size 2^ and "rotate" all links 
going from sites x G x' in a fixed r-direction with the same angle O^': 

U,,r~^Kr = 9:Rx9xU,,r, (82) 

with 

Rx{nx',Ox') = coslOx'/'i) + ism{6x'/2) fix'- (y . (83) 

The gauge transformation g is obtained by imposing the Coulomb gauge condition 
on slices A[ as defined above. We now consider the special case where the directions 
fix' of the embedded ?7(l)-subgroups are independent of the block x', i.e. fix' = n 
for all x' G Ai. Then we get a conditional Hamiltonian 'Hi{9) by substituting 
the "rotated" gauge field U' in the fundamental Hamiltonian. By iterating this 
procedure one gets conditional Hamiltonians Ti.k{0) on coarser layers A^,. The point 
is that in the special case considered here T-Ck{6) always stays of the form 

-'HkiO) = ^ J2 {Px'Cos^iOx')+P:icosiex')smiex') + P'Jsm^Ox')} 

+ ^ E E {(^x',^. <^os{ex') cosiOx'+f,) + 1311^^ cos(^,0 sin(^,,+^) 

+ (311^ sin{ex') cos{ex'+f,) + sin(^^0 sin(^,.+^) } 

+ const , (84) 

with space dependent couplings that can be recursively calculated from the couplings 
defined on the next finer layer Afc_i. Therefore a W-cycle is possible. Of course, 
this implementation is also possible with three dimensional blocks. 
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8 Summary 



We have presented a simple yet accurate formula that expresses acceptance rates 
for nonlocal update algorithms in terms of one single parameter (or two in the case 
of nonabelian gauge theory) entering the quantity hi. This parameter is easy to 
compute, e.g. by Monte Carlo simulations on a small lattice. We encountered two 
classes of models. For Sine Gordon, 0^ theory and SU{2) lattice gauge theory, s 
had to be rescaled like for piecewise constant and for smooth kernels, whereas 
for massless free field theory, the XY model, the 0{N) nonlinear cr-model and U{1) 
lattice gauge theory, one can achieve Ls-independent acceptance rates by choosing 
smooth kernels. 

We can compare the behavior of the acceptance rates in interacting models with 
free field theory, where CSD is known to be eliminated by a multigrid algorithm. In 

order to do this we presented a study of the infiuence of the coarse-to-fine interpo- 
lation on the acceptance rates in free field theory. 

The results of the comparison arc consistent with the following rule: For an 
interacting model, sufficiently high acceptance rates for a complete elimination of 
CSD can only be expected if hi — -|- st/j) — 'H{4>)) contains no algorithmic 

"mass" term ~ With the help of this rule it is possible to decide whether 

a given statistical model is a natural candidate for multigrid Monte Carlo or not. 

The kinematical mechanism that leads to a failure of multigrid algorithms is well 
described by our analysis. Wc hope that a better understanding can lead to improved 
multigrid algorithms that can overcome kinematical obstructions stemming from an 
algorithmic "mass" term. 

The acceptance rates of our proposal for nonlocal updates in SU{2) lattice gauge 
theory were investigated in detail. Here we found that an algorithmic "mass" term 
generated by the disorder in the gauge field suppresses the acceptance rates on large 
blocks. From this study we do not expect that our algorithm will have a chance to 
overcome CSD. However, we believe that compared to local Metropolis algorithms 
there will be an acceleration of the Monte Carlo dynamics by a constant factor 
(depending on the details of the implementation). We think that the best method 
for practical purposes would be an explicit multigrid implementation using piecewise 
constant kernels and a W-cycle. An implementation and test is planned. The crucial 
question is, of course, whether one is able to beat the overrelaxation algorithm. 
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Figure Captions 



FIG.l: Metropolis step sizes e{LB) for massless free field theory in two dimensions, 
512^-lattice. e{LB) is choosen in such a way that always Pace = 0.5 holds. Symbols: 
full circles: ip^^"''^^ ^ full triangles: empty circles: empty triangles: ■j/'**"'^, 

full squares: -^''^^^"''^ empty squares: -0'^°"**. Lines are only drawn to guide the eye. 

FIG. 2: Metropolis step sizes e{LB) for massless free field theory in four dimensions, 
64^-lattice. s{Lb) is choosen in such a way that always Pace — 0.5 holds. Symbols: 
full circles: ij)^^'^'^'^^ full triangles: -0*™"'=. empty circles: -0™", empty triangles: ■0*'"^, 
full squares: ip''^^'^^'^ ^ empty squares: ■0'=°"^*. Lines are only drawn to guide the eye. 

FIG.3: Q{s) for piecewise constant kernels in the 2-dimensional Sine Gordon model, 
p = 39.478, C = 1- From top to bottom: Lb = 4,8, 16,32 on a 16^, 32^,642, 128^ 
lattice, respectively. Points with error bars: Monte Carlo results, hues: analytical 
results. 

FIG. 4: Q{s) for piecewise constant kernels in the 2-dimensional XY model, /3 = 1.2. 
From top to bottom: Lb = 4,8,16 on a 16^,32^,64^ lattice, respectively. Points 
with error bars: Monte Carlo results, lines: analytical results. 

FIG. 5: Q{s) for piecewise constant kernels in the 2-dimensional 0^ theory, = 
-0.56, Xo = 2.4. From top to bottom: Lb = 4, 8, 16 on a 16^ 32^, 64^ lattice, 
respectively. Points with error bars: Monte Carlo results, lines: analytical results. 

FIG. 6: fl{s) in 4-dimensional SU{2) lattice gauge theory using piecewise constant 
kernels, j3 = 2.6 on a 20^-lattice. From top to bottom: Lb = 2,4,8,16. Points 
with error bars: Monte Carlo results, lines: analytical results using rriD and P from 
Monte Carlo (errors smaller than line width). 

FIG. 7: Comparison of disorder and surface effects for 4-dimensional SU{2) lattice 
gauge theory using piecewise constant kernels on an 8^-block, P = 2.6 on a 20^- 
lattice. Sohd line: (disorder effects), dashed line: /ii,p(s) (surface effects). 

FIG. 8: Comparison of disorder and surface effects for the 4-dimensional SU{2) 
lattice gauge theory using smooth -0**"'^ kernels on a 8'^-block, (3 = 2.6 on a 20^- 
lattice, quadratic approximation used. Solid line: /ii,^(s) (disorder effects), dashed 
line: /ii^p_^(s) (surface effects). 
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